#ifndef CURVE_FITTING_H
#define CURVE_FITTING_H

#include <vector>
#include <cmath>
#include <iostream>
#include <sstream>
#include "Function.hpp"
#include "matplotlibcpp.h"
#include "spline_functions.h"

using namespace std;
namespace plt = matplotlibcpp;

// 先写2维的标准t参数
void Curve_fitting_2_t( const Function&f, int N, double t_min, double t_max, int type = 1, double mx1 = 0, double mxN = 0, double my1 = 0, double myN = 0) {
    // f是用参数表示的曲线，N是给定分段数，剩下是需要拟合的曲线的范围，type标记拟合时的边界条件（周期为0，完全为1）

    double gap = static_cast<double>(t_max - t_min)/N;
    std::vector<double> knots_x, knots_y;
    std::vector<double> coefficients_x, coefficients_y;
    double t = t_min;
    for (int i = 0; i < N+1; i++ ) {
        if( i==N ){
            t = t_max;
        }
        dot point = f(t);
        knots_x.push_back(t);
        knots_y.push_back(t);
        coefficients_x.push_back(point.x);
        coefficients_y.push_back(point.y);
        t += gap;        
    }
    BSpline* Spline_B_x = nullptr;
    BSpline* Spline_B_y = nullptr;

    if(type == 0){
        Spline_B_x = new PeriodicCubicSpline_B(3, knots_x, coefficients_x);
        Spline_B_y = new PeriodicCubicSpline_B(3, knots_y, coefficients_y);
    }
    else{
        Spline_B_x = new ClampedCubicSpline_B(3, knots_x, coefficients_x, mx1, mxN);
        Spline_B_y = new ClampedCubicSpline_B(3, knots_y, coefficients_y, my1, myN);
    }

    vector<double> x_values, y_values;
    int M = 100;
    gap = static_cast<double>(t_max - t_min)/M;
    t = t_min;
    for (int i = 0; i < M; i++ ) { // 为什么不能是M，因为会积累误差导致无法正好等于t_max，所以在后面单独添加最后一个点
        x_values.push_back(Spline_B_x->evaluate(t));
        y_values.push_back(Spline_B_y->evaluate(t));
        t += gap;        
    } 
    // x_values.push_back(Spline_B_x->evaluate(t_max));
    // y_values.push_back(Spline_B_y->evaluate(t_max)); 
          
    plt::named_plot("Spline_Curve_2_t", x_values, y_values, "r-");
    plt::xlabel("X");
    plt::ylabel("Y");
    plt::grid(true);
    plt::legend();
    //plt::show();
}
// 2维的弦长参数
void Curve_fitting_2_s( const Function&f, int N, double t_min, double t_max, int type = 1, double mx1 = 0, double mxN = 0, double my1 = 0, double myN = 0) { 
    // f是用参数表示的曲线，N是给定分段数，剩下是需要拟合的曲线的范围，type标记拟合时的边界条件（周期为0，完全为1）

    double gap = static_cast<double>(t_max - t_min)/N;
    std::vector<double> knots_x, knots_y;
    std::vector<double> coefficients_x, coefficients_y;
    double t = t_min;
    double s = 0;
    dot point_1 = f(t);
    dot point_2;
    for (int i = 0; i < N+1; i++ ) {
        if( i==N ){
            t = t_max;
        }
        point_2 = f(t);
        double delta_s = pow(pow((point_2.x-point_1.x),2)+pow((point_2.y-point_1.y),2),0.5);
        s += delta_s;
        knots_x.push_back(s);
        knots_y.push_back(s);
        coefficients_x.push_back(point_2.x);
        coefficients_y.push_back(point_2.y);
        point_1 = point_2;
        t += gap;        
    }
    double s_max = s;
    BSpline* Spline_B_x = nullptr;
    BSpline* Spline_B_y = nullptr;

    if(type == 0){
        Spline_B_x = new PeriodicCubicSpline_B(3, knots_x, coefficients_x);
        Spline_B_y = new PeriodicCubicSpline_B(3, knots_y, coefficients_y);
    }
    else{
        Spline_B_x = new ClampedCubicSpline_B(3, knots_x, coefficients_x, mx1, mxN);
        Spline_B_y = new ClampedCubicSpline_B(3, knots_y, coefficients_y, my1, myN);
    }

    vector<double> x_values, y_values;
    int M = 100;
    gap = static_cast<double>(s_max)/M;
    s = 0;
    for (int i = 0; i < M; i++ ) { // 为什么不能是M，因为会积累误差导致无法正好等于t_max，所以在后面单独添加最后一个点
        x_values.push_back(Spline_B_x->evaluate(s));
        y_values.push_back(Spline_B_y->evaluate(s));
        s += gap;        
    } 
    // x_values.push_back(Spline_B_x->evaluate(s_max));
    // y_values.push_back(Spline_B_y->evaluate(s_max)); 
          
    plt::named_plot("Spline_Curve_2_s", x_values, y_values, "g-");
    plt::xlabel("X");
    plt::ylabel("Y");
    plt::grid(true);
    plt::legend();
    //plt::show();
}

// 球面的t参数
void Curve_fitting_Sphere_t(double R, const Function&f, int N, double t_min, double t_max, int type = 1, double mx1 = 0, double mxN = 0, double my1 = 0, double myN = 0) {
    // f是用参数表示的曲线，N是给定分段数，剩下是需要拟合的曲线的范围，type标记拟合时的边界条件（周期为0，完全为1）, R 为球面半径
    // 以下的x,y实际上是纬度v和经度u
    double gap = static_cast<double>(t_max - t_min)/N;
    std::vector<double> knots_x, knots_y;
    std::vector<double> coefficients_x, coefficients_y;
    double t = t_min;
    for (int i = 0; i < N+1; i++ ) {
        if( i==N ){
            t = t_max;
        }
        dot point = f(t);
        knots_x.push_back(t);
        knots_y.push_back(t);
        coefficients_x.push_back(point.x);
        coefficients_y.push_back(point.y);
        t += gap;        
    }
    BSpline* Spline_B_x = nullptr;
    BSpline* Spline_B_y = nullptr;

    if(type == 0){
        Spline_B_x = new PeriodicCubicSpline_B(3, knots_x, coefficients_x);
        Spline_B_y = new PeriodicCubicSpline_B(3, knots_y, coefficients_y);
    }
    else{
        Spline_B_x = new ClampedCubicSpline_B(3, knots_x, coefficients_x, mx1, mxN);
        Spline_B_y = new ClampedCubicSpline_B(3, knots_y, coefficients_y, my1, myN);
    }

    vector<double> x_values, y_values, z_values;
    int M = 100;
    gap = static_cast<double>(t_max - t_min)/M;
    t = t_min;
    for (int i = 0; i < M; i++ ) { // 为什么不能是M，因为会积累误差导致无法正好等于t_max，所以在后面单独添加最后一个点
        x_values.push_back(R * sin(Spline_B_x->evaluate(t)) * cos(Spline_B_y->evaluate(t)));
        y_values.push_back(R * sin(Spline_B_x->evaluate(t)) * sin(Spline_B_y->evaluate(t)));
        z_values.push_back(R * cos(Spline_B_x->evaluate(t)));
        t += gap;        
    } 
    x_values.push_back(R * sin(Spline_B_x->evaluate(t_max)) * cos(Spline_B_y->evaluate(t_max)));
    y_values.push_back(R * sin(Spline_B_x->evaluate(t_max)) * sin(Spline_B_y->evaluate(t_max)));
    z_values.push_back(R * cos(Spline_B_x->evaluate(t_max)));
          
    plt::plot3(x_values, y_values, z_values); 
    plt::xlabel("X");
    plt::ylabel("Y");
    plt::grid(true);
    //plt::legend();
    plt::show();
}
// 球面的弦长参数
void Curve_fitting_Sphere_s(double R, const Function&f, int N, double t_min, double t_max, int type = 1, double mx1 = 0, double mxN = 0, double my1 = 0, double myN = 0) { 
    // f是用参数表示的曲线，N是给定分段数，剩下是需要拟合的曲线的范围，type标记拟合时的边界条件（周期为0，完全为1）
    // 以下的x,y实际上是纬度v和经度u
    double gap = static_cast<double>(t_max - t_min)/N;
    std::vector<double> knots_x, knots_y;
    std::vector<double> coefficients_x, coefficients_y;
    double t = t_min;
    double s = 0;
    dot point_1 = f(t);
    dot point_2;
    for (int i = 0; i < N+1; i++ ) {
        if( i==N ){
            t = t_max;
        }
        point_2 = f(t);
        double delta_s = pow(pow((point_2.x-point_1.x),2)+pow((point_2.y-point_1.y),2),0.5);
        s += delta_s;
        knots_x.push_back(s);
        knots_y.push_back(s);
        coefficients_x.push_back(point_2.x);
        coefficients_y.push_back(point_2.y);
        point_1 = point_2;
        t += gap;        
    }
    double s_max = s;
    BSpline* Spline_B_x = nullptr;
    BSpline* Spline_B_y = nullptr;

    if(type == 0){
        Spline_B_x = new PeriodicCubicSpline_B(3, knots_x, coefficients_x);
        Spline_B_y = new PeriodicCubicSpline_B(3, knots_y, coefficients_y);
    }
    else{
        Spline_B_x = new ClampedCubicSpline_B(3, knots_x, coefficients_x, mx1, mxN);
        Spline_B_y = new ClampedCubicSpline_B(3, knots_y, coefficients_y, my1, myN);
    }

    vector<double> x_values, y_values, z_values;
    int M = 100;
    gap = static_cast<double>(s_max)/M;
    s = 0;
    for (int i = 0; i < M; i++ ) { // 为什么不能是M，因为会积累误差导致无法正好等于t_max，所以在后面单独添加最后一个点
        x_values.push_back(R * sin(Spline_B_x->evaluate(s)) * cos(Spline_B_y->evaluate(s)));
        y_values.push_back(R * sin(Spline_B_x->evaluate(s)) * sin(Spline_B_y->evaluate(s)));
        z_values.push_back(R * cos(Spline_B_x->evaluate(s)));
        s += gap;        
    } 
    x_values.push_back(R * sin(Spline_B_x->evaluate(s_max)) * cos(Spline_B_y->evaluate(s_max)));
    y_values.push_back(R * sin(Spline_B_x->evaluate(s_max)) * sin(Spline_B_y->evaluate(s_max)));
    z_values.push_back(R * cos(Spline_B_x->evaluate(s_max)));
          
    plt::plot3(x_values, y_values, z_values); 
    plt::xlabel("X");
    plt::ylabel("Y");
    plt::grid(true);
    //plt::legend();
    plt::show();
}

// 绘制真实曲线
void plot_exact_2 (const Function& f, const double& t1, const double& t2){
    std::vector<double> x_exact, y_exact;
    int N = 1000;
    double h = (t2 - t1) /N;
    for (int i = 0; i <= N; ++i) {
        x_exact.push_back(f(t1+i*h).x);
        y_exact.push_back(f(t1+i*h).y);
    }
    // 使用 matplotlib-cpp 绘图
    plt::named_plot("Exact Curve", x_exact, y_exact, "b--");
    plt::legend();
    plt::show();
}
void plot_exact_3 (const Function_3& f, const double& t1, const double& t2){
    std::vector<double> x_exact, y_exact, z_exact;
    int N = 1000;
    double h = (t2 - t1) /N;
    for (int i = 0; i <= N; ++i) {
        x_exact.push_back(f(t1+i*h).x);
        y_exact.push_back(f(t1+i*h).y);
        z_exact.push_back(f(t1+i*h).z);
    }

    // 使用 matplotlib-cpp 绘图
    plt::plot3(x_exact, y_exact, z_exact); 
    plt::legend();
    //plt::show();
}

/*
// 3维的标准t参数
void Curve_fitting_3_t( const Function_3&f, int N, double t_min, double t_max) { 
    // f是用参数表示的曲线，N是给定分段数，剩下是需要拟合的曲线的范围

    double gap = static_cast<double>(t_max - t_min)/N;
    std::vector<double> knots_x, knots_y, knots_z;
    std::vector<double> coefficients_x, coefficients_y, coefficients_z;
    double t = t_min;
    for (int i = 0; i < N+1; i++ ) {
        if( i==N ){
            t = t_max;
        }
        dot_3 point = f(t);
        knots_x.push_back(t);
        knots_y.push_back(t);
        knots_z.push_back(t);
        coefficients_x.push_back(point.x);
        coefficients_y.push_back(point.y);
        coefficients_z.push_back(point.z);
        t += gap;        
    }
    PeriodicCubicSpline_B periodicSpline_B_x(3, knots_x, coefficients_x);
    PeriodicCubicSpline_B periodicSpline_B_y(3, knots_y, coefficients_y);
    PeriodicCubicSpline_B periodicSpline_B_z(3, knots_z, coefficients_z);

    vector<double> x_values, y_values, z_values;
    // vector<double> x_exacts, y_exacts, z_exacts;
    int M = 100;
    gap = static_cast<double>(t_max - t_min)/M;
    t = t_min;
    for (int i = 0; i < M; i++ ) { // 为什么不能是M，因为会积累误差导致无法正好等于t_max，所以在后面单独添加最后一个点
        x_values.push_back(periodicSpline_B_x.evaluate(t));
        y_values.push_back(periodicSpline_B_y.evaluate(t));
        z_values.push_back(periodicSpline_B_z.evaluate(t));
        // x_exacts.push_back(f(t).x);
        // y_exacts.push_back(f(t).y);
        // z_exacts.push_back(f(t).z);
        t += gap;        
    } 
    x_values.push_back(periodicSpline_B_x.evaluate(t_max));
    y_values.push_back(periodicSpline_B_y.evaluate(t_max)); 
    z_values.push_back(periodicSpline_B_z.evaluate(t_max));
    // x_exacts.push_back(f(t_max).x);
    // y_exacts.push_back(f(t_max).x);
    // z_exacts.push_back(f(t_max).x);
          
    // 使用plot3绘制三维曲线
    plt::plot3(x_values, y_values, z_values); 
    // plt::plot3(x_exacts, y_exacts, z_exacts);
    plt::xlabel("X");
    plt::ylabel("Y");
    // plt::zlabel("Z");
    plt::grid(true);
    plt::show();
}
// 3维的弦长参数
void Curve_fitting_3_s( const Function_3&f, int N, double t_min, double t_max) { 
    // f是用参数表示的曲线，N是给定分段数，剩下是需要拟合的曲线的范围，

    double gap = static_cast<double>(t_max - t_min)/N;
    std::vector<double> knots_x, knots_y, knots_z;
    std::vector<double> coefficients_x, coefficients_y, coefficients_z;
    double t = t_min;
    double s = 0;
    dot_3 point_1 = f(t);
    dot_3 point_2;
    for (int i = 0; i < N+1; i++ ) {
        if( i==N ){
            t = t_max;
        }
        point_2 = f(t);
        double delta_s = pow(pow((point_2.x-point_1.x),2)+pow((point_2.y-point_1.y),2)+pow((point_2.z-point_1.z),2),0.5);
        s += delta_s;
        knots_x.push_back(s);
        knots_y.push_back(s);
        knots_z.push_back(s);
        coefficients_x.push_back(point_2.x);
        coefficients_y.push_back(point_2.y);
        coefficients_z.push_back(point_2.z);
        point_1 = point_2;
        t += gap;        
    }
    double s_max = s;
    PeriodicCubicSpline_B periodicSpline_B_x(3, knots_x, coefficients_x);
    PeriodicCubicSpline_B periodicSpline_B_y(3, knots_y, coefficients_y);
    PeriodicCubicSpline_B periodicSpline_B_z(3, knots_z, coefficients_z);

    vector<double> x_values, y_values, z_values;
    int M = 100;
    gap = static_cast<double>(s_max)/M;
    s = 0;
    for (int i = 0; i < M; i++ ) { // 为什么不能是M，因为会积累误差导致无法正好等于t_max，所以在后面单独添加最后一个点
        x_values.push_back(periodicSpline_B_x.evaluate(s));
        y_values.push_back(periodicSpline_B_y.evaluate(s));
        z_values.push_back(periodicSpline_B_z.evaluate(s));
        s += gap;        
    } 
    x_values.push_back(periodicSpline_B_x.evaluate(s_max));
    y_values.push_back(periodicSpline_B_y.evaluate(s_max)); 
    z_values.push_back(periodicSpline_B_z.evaluate(s_max));

    plt::plot3(x_values, y_values, z_values);
    plt::xlabel("X");
    plt::ylabel("Y");
    plt::grid(true);
    plt::show();
}
*/

#endif // CURVE_FITTING_H